7.1.1 R包地图底图

## 获取shp地图：
newmap <- getMap(resolution = "coarse")  
plot(newmap)
data(coastsCoarse) 
data(countriesLow) 

## 直接绘图，然后添加shp边界：
hh = read.csv("./dhfuh.csv")
plot(hh)
maps::maps(add= T)

## 添加box获取地图边界：
data(wrld_simpl)

# Plot the base map
plot(wrld_simpl, 
     xlim = c(min.lon, max.lon),
     ylim = c(min.lat, max.lat),
     axes = TRUE, 
     col = "grey95")、
box()

7.12 leaflet

library(wallace)
leaflet()%>%setView(lng=116.38,lat=39.9,zoom=3)%>%
  addTiles()%>%addProviderTiles("Esri.WorldImagery")

## 查看范围：
data(wrld_simpl)

## 使用mao(add =T)进行快速可视化：
all <- read.csv("D:/xh2/f1_points/xhas.csv")
plot(all$longitude,all$latitude)
map(add=T)

## ggplot代码的快速并列可视化：
grid.arrange(plot1, plot2, ncol = 2)

​```{r, fig.width=8}
plot1 = ggplot(combined_dat, aes(y = ..density..)) + geom_histogram(aes(x = Sepal.Width, fill = Type), bins = 20) +
  facet_wrap(~Type) +
  ggtitle("Distribution of Sepal Width", "Biased resample vs Original sample")
plot2 = ggplot(combined_dat, aes(y = ..density..)) + geom_histogram(aes(x = Sepal.Length, fill = Type), bins = 20) +
  facet_wrap(~Type) +
  ggtitle("Distribution of Sepal Length", "Biased resample vs Original sample")
grid.arrange(plot1, plot2, ncol = 2)
​

案例：

加载R包：

library(raster) library(dismo) library(tidyverse) library(SDMtune) library(rgdal)

设置相对路径：

setwd("E:/sjdata/")

加载物种分布数据：

待比较中四个种间数据：

加载西藏沙棘：

hthi <- read.csv("./F1_SP/Hthi3.csv")

加载江孜沙棘：

hgya <- read.csv("./F1_SP/Hgya3.csv")

加载柳叶沙棘：

hsal <- read.csv("F1_SP/Hsal3.csv")

加载肋果沙棘：

hsneu <- read.csv("./F1_SP/Hsneu3.csv") hsste <- read.csv("./F1_SP/Hsste3.csv") hse <- rbind(hsneu,hsste ) hse2 <- rep("Hse",dim(hse)[1]) hse <- cbind(hse2,hse[,2:3]) names(hse)[1] <- "name"

汇总分布数据：

alls <- rbind(hthi,hgya,hsal,hse) summary(all) head(hthi) head(hgya) head(hse)

加载鼠李沙棘的物种分布数据：

中亚沙棘 stur

中国沙棘 ssin

高加索沙棘 scau

云南沙棘 syun

蒙古沙棘 smon

卧龙沙棘 swol

海滨沙棘 srha

克尔巴阡山沙棘 scar

溪生沙棘 sflu

stur <- read.csv("./F1_SP/stur3.csv") ssin <- read.csv("./F1_SP/ssin3.csv") scau <- read.csv("./F1_SP/scau3.csv") syun <- read.csv("./F1_SP/syun3.csv") smon <- read.csv("./F1_SP/smon3.csv") srha <- read.csv("./F1_SP/srha3.csv") scar <- read.csv("./F1_SP/scar3.csv") sflu <- read.csv("./F1_SP/sflu5.csv")

all <- rbind(hthi,hgya,hse) summary(all)

long :Min: 77.27 Max:103.46

lat :min:27.72 max:39.91

构建5度buffer：

frange <- function(x){ occs <- x xmin <- min(occs$longitude) xmax <- max(occs$longitude) ymin <- min(occs$latitude) ymax <- max(occs$latitude) bb <- matrix(c(xmin-4, xmin-4, xmax+4, xmax+4, xmin-4, ymin-3, ymax+4, ymax+4, ymin-3, ymin-3), ncol=2) bgExt <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(bb)), 1))) return(bgExt) }

rang <- frange(all) rang

导出数据分布范围：

df <- data.frame(ID=character(), stringsAsFactors=FALSE ) for (i in rang@polygons ) { df <- rbind(df, data.frame(ID=i@ID, stringsAsFactors=FALSE)) } row.names(df) <- df$ID rang <- SpatialPolygonsDataFrame(rang, df)

可视化种间分布点，使用leaflet：

加载属：

all3 <- rbind(hthi,hgya,hsal,hse) %>% data.frame(.) head(all3) attach(all3) library(leaflet)

iconList = awesomeIconList( "Hthi" = makeAwesomeIcon(icon ='ios-close',markerColor = "blue"), "Hgya" = makeAwesomeIcon(icon ='ios-close',markerColor = "green"), "Hsal" = makeAwesomeIcon(icon ='ios-close',markerColor = "red"), "Hse" = makeAwesomeIcon(icon ='ios-close',markerColor = "black")) pal <- colorNumeric( palette = "YlGnBu", domain = countries$gdp_md_est )

all3$name <- factor(all3$name)

yyIcons <- iconList( 中国医科院所属医院 = makeIcon(“https://img-blog.csdn.net/20161015181859390”, iconWidth =32, iconHeight = 32), 北京区县属医院 = makeIcon(“https://img-blog.csdn.net/20161015182217958”, iconWidth =32, iconHeight = 32),

leaflet(bj3H) %>%addProviderTiles(“CartoDB.Positron”) %>% addMarkers(icon = ~yyIcons[fl],popup=~fl)

iconList = iconList( Hthi = makeIcon(icon ='ios-close',markerColor = "blue"), Hgya = makeIcon(icon ='ios-close',markerColor = "green"), Hsal = makeIcon(icon ='ios-close',markerColor = "red"), Hse = makeIcon(icon ='ios-close',markerColor = "black"))

有待改进的地方：

参考url：https://www.giserdqy.com/secdev/leaflet/25986/

https://www.baidu.com/s?ie=UTF-8&wd=R%E8%AF%AD%E8%A8%80%E5%9C%A8%E7%BA%BF%E5%9C%B0%E5%9B%BE%E7%A5%9E%E5%99%A8%EF%BC%9ALeaflet%20for%20R%E5%8C%85

添加图例：

p = colorFactor(palette = c("white","green","blue"),domain = c("Office","college","home"),ordered = T) my_adds <- leaflet() %>% addTiles() %>% addCircleMarkers(lng=75.6944, lat=11.9811,popup = "My Office",color = "white") %>% addCircleMarkers(lng=77.7141, lat=12.9675,popup = "My college",color = "green") %>% addCircleMarkers(lng=77.669907, lat=13.0118285,popup = "My Home") %>% addLegend(position = "topright",pal = p, values = c("office","home","college"),title = "Address") my_adds

国界分布简图：Esri.WorldStreetMap

leaflet(data = all3[,-1]) %>% addTiles() %>% addMarkers(~longitude, ~latitude)%>% addAwesomeMarkers(lng=~longitude,lat=~latitude,icon = ~iconList[name], label= ~as.character(name)) %>% addProviderTiles("Esri.WorldStreetMap")

换卫星图：

leaflet(data = all3[,-1]) %>% addTiles() %>% addMarkers(~longitude, ~latitude)%>% addAwesomeMarkers(lng=~longitude,lat=~latitude,icon = ~iconList[name], label= ~as.character(name)) %>% addProviderTiles("Stamen.Terrain") %>% addPolygons(data =rang,fill = F)

pal <- colorFactor(c("blue","green","red","black"), all3, levels = NULL, ordered = FALSE, na.color = "#808080", alpha = FALSE, reverse = FALSE)

加载一个研究范围：

pal <- colorNumeric( palette = "YlGnBu", domain = countries$gdp_md_est ) map %>% addPolygons(stroke = FALSE, smoothFactor = 0.2, fillOpacity = 1, color = ~pal(gdp_md_est) ) %>% addLegend("bottomright", pal = pal, values = ~gdp_md_est, title = "Est. GDP (2010)", labFormat = labelFormat(prefix = "$"), opacity = 1 )

各种地图底图

# OpenStreetMap.Mapnik

OpenStreetMap.BlackAndWhite

OpenStreetMap.DE

OpenStreetMap.France

OpenStreetMap.HOT

OpenTopoMap

Thunderforest.OpenCycleMap

Thunderforest.Transport

Thunderforest.TransportDark

Thunderforest.SpinalMap

Thunderforest.Landscape

Thunderforest.Outdoors

Thunderforest.Pioneer

OpenMapSurfer.Roads

OpenMapSurfer.Grayscale

Hydda.Full

Stamen.Toner

Stamen.TonerBackground

Stamen.TonerLite

Stamen.Watercolor

Stamen.Terrain

Stamen.TerrainBackground

Stamen.TopOSMRelief

Esri.WorldStreetMap

Esri.DeLorme

Esri.WorldTopoMap

Esri.WorldImagery

Esri.WorldTerrain

Esri.WorldShadedRelief

Esri.WorldPhysical

Esri.OceanBasemap

Esri.NatGeoWorldMap

Esri.WorldGrayCanvas

MtbMap

CartoDB.Positron

CartoDB.PositronNoLabels

CartoDB.PositronOnlyLabels

CartoDB.DarkMatter

CartoDB.DarkMatterNoLabels

CartoDB.DarkMatterOnlyLabels

HikeBike.HikeBike

HikeBike.HillShading

NASAGIBS.ModisTerraTrueColorCR

NASAGIBS.ModisTerraBands367CR

NASAGIBS.ViirsEarthAtNight2012

NASAGIBS.ModisTerraLSTDay

NASAGIBS.ModisTerraSnowCover

NASAGIBS.ModisTerraAOD

NASAGIBS.ModisTerraChlorophyll

### 优秀的leaflet案例

```r
### 在浏览器中打开：
## 这种leaflet的使用方法值得学习：
m <- leaflet()
m <- addTiles(m)
m <- addProviderTiles(m, "Esri.OceanBasemap")

cols <- c("red", "navy")
m <- addCircleMarkers(m, 
                      lng = hurricanes$FirstLon, 
                      lat = hurricanes$FirstLat,
                      radius = 2.5,
                      color = cols[hurricanes$Type+1],
                      popup = paste('Year:', as.character(hurricanes$Year)))
m <- addLegend(m,
               "topright", 
               colors = cols, 
               labels = c('tropical', 'non-tropical'),
               title = "Type of Hurricane",
               opacity = 1)
m
#

7.1.3 卫星地图可视化：

## 卫星地图可视化：
e = extent(chuck.sp)
buf = .5
chuck.df<- as.data.frame(chuck.sp)

library(ggmap)
myMap <- get_stamenmap(bbox = c(left = e[1]-buf,
                                bottom = e[3] -buf,
                                right = e[2]+buf,
                                top = e[4] +buf),
                       maptype = "terrain",
                       crop = FALSE,
                       zoom = 6)
                       # plot map
ggmap(myMap) +  geom_point(aes(x = lon, y = lat), data = chuck.df, alpha = .5)


## 谷歌地图可视化：
library(RgoogleMaps)
lat <- c(48,64) #define our map's ylim
lon <- c(-140,-110) #define our map's xlim
center = c(mean(lat), mean(lon))  #tell what point to center on
zoom <- 5  #zoom: 1 = furthest out (entire globe), larger numbers = closer in
terrmap <- GetMap(center=center, zoom=zoom, maptype= "terrain", destfile = "terrain.png") #lots of visual options, just like google maps: maptype = c("roadmap", "mobile", "satellite", "terrain", "hybrid", "mapmaker-roadmap", "mapmaker-hybrid")

7.1.4 卫星地图可视化：

## 使用mapview进行可视化
###### f2:分布图可视化：####
t1 <- occ[,2:3] %>% data.frame()
library(mapview)
# install.packages("mapview")
library(sp)
coordinates(t1 ) <- ~long+lat
proj4string(t1) <- CRS("+init=epsg:4326")
# plot data
mapView(t1)